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ABSTRACT 

We analyzed the soft X-ray light curves from the Geostationary Operational Environmental 
Satellites (GOES) over the last 37 years (1975-2011) and measured with an automated flare 
detection algorithm over 300,000 solar flare events (amounting to w 5 times higher sensitivity 
than the NOAA flare catalog). We find a powerlaw slope of otp = 1.98±0.11 for the (background- 
subtracted) soft X-ray peak fluxes that is invariant through three solar cycles and agrees with 
the theoretical prediction ap = 2.0 of the fractal-diffusive self-organized criticality (FD-SOC) 
model. For the soft X-ray flare rise times we find a powerlaw slope of cut = 2.02 ± 0.04 during 
solar cycle minima years, which is also consistent with the prediction cut = 2.0 of the FD- 
SOC model. During solar cycle maxima years, the powerlaw slope is steeper in the range of 
ctr ~ 2.0 — 5.0, which can be modeled by a solar cycle- dependent flare pile-up bias effect. These 
results corroborate the FD-SOC model, which predicts a powerlaw slope of qe = 1.5 for flare 
energies and thus rules out significant nanoflare heating. While the FD-SOC model predicts the 
probability distribution functions of spatio-temporal scaling laws of nonlinear energy dissipation 
processes, additional physical models are needed to derive the scaling laws between the geometric 
SOC parameters and the observed emissivity in different wavelength regimes, as we derive here 
for soft X-ray emission. The FD-SOC model yields also statistical probabilities for solar flare 
forecasting. 

Subject headings: Sun: corona — Sun: X-rays — Sun: flares — methods: statistical 



1. INTRODUCTION 

There are a number of intriguing questions that have been raised about statistics of solar flares: (1) Why 
do they show the ubiquitous powerlaws in their occurrence frequency distributions? (2) Do we understand 
the numerical values of the powerlaw slopes? (3) Do the powerlaw slopes vary with the solar cycle? (4) Is 
the solar corona heated by nanoflares ? (5) Can the flare statistics be explained in terms of the self-organized 
criticality (SOC) concept? (6) What are the consequences of SOC models on our physical understanding of 
solar flares? (7) How can we improve statistical solar flare prediction? (8) What are the largest expected 
solar flares in history and future? In this paper we address these questions by conducting statistics of solar 
flares using the longest uniform dataset we have available, namely the soft X-ray light curves from the 
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Geostationary Operational Environmental Satellites (GOES), which cover over 37 years during the period 
of 1974-2012. 

The GOES program consists of a series of geostationary satellites (orbiting the Earth at a height of 
35,790 km), which overlap in time so that there are always one to three spacecraft present and guarantee an 
essentially uninterrupted time series of solar soft X-ray fluxes, besides continuous meteorological observations 
of the Earth. The GOES-1 satellite was launched on October 16, 1974, and GOES-2 and GOES-3 followed in 
1977 and 1978. In the meantime the series continued up to GOES-15, launched on March 4, 2010, while the 
future satellites GOES-R and GOES-S (with soft X-ray imaging capabilities) are in the queue for a launch in 
2015 and 2017, respectively. Operational and technical details of GOES satellites can be gleaned from Grubb 
(1975), Donnelly et al. (1977), Bouwer et al. (1982), Thomas et al. (1985), Kahler and Kreplin (1991), Garcia 
(1994), Hill (2005), Pizzo et al. (2005), Shing et al. (1999), White et al. (2005), Neupert (2011), the National 
Oceanographic and Atmospheric Administration (NOAA) website http://www.oso.noaa.gov/goes/, or the 



NASA website http://goespoes.gsfc.nasa.gov/project/index.html For our study we are concerned with the 



soft X-ray light curves, which are recorded in two energy channels, i.e., in the softer energy range of 1-8 A 
and in the harder energy range of 0.5-4 A. Both light curves are available with a cadence of 3 s, and with a 
cadence of 2 s after 2009 Dec 1. 

Solar flare statistics can be generated by detection of sudden impulsive increases of the soft X-ray 
intensity in light curves. NOAA publishes solar flare catalogs that are issued on a daily basis. Here we 
develop an automated flare detection algorithm that allows us to analyze the entire 37-year time series in an 
objective way without human subjectivity. We are interested in the mathematical function of the probability 
distributions of solar flare parameters. It has been established earlier that the size distribution of soft X-ray 
peak fluxes of solar flares have a powerlaw-like form, with a powerlaw slope in the range of ap fts 1.6 — 2.1 
(Hudson et al. 1969, Drake et al. 1971, Shimizu 1995, Lee et al. 1995, Feldman et al. 1997, Shimojo and 
Shibata 1999, Veronig et al. 2002a, b, Yashiro et al. 2006). In this study we will measure this parameter 
with the largest available statistics, including over 300,000 detected solar flare events, and will study its 
variation during three solar cycles. Further we will investigate measurement uncertainties of previously 
analyzed distributions and compare them with the results from this study in order to establish the most 
accurate value. 

The appearance of powerlaw-like size distribution functions has generally been linked to the theory of 
self-organized criticality (SOC) (or alternatively to turbulence), which predicts scale- free powerlaw distri- 
butions (Bak et al. 1987, 1988). A prominent paradigm of the SOC concept is a sandpile that produces 
a scale-free distribution of avalanche sizes, once it reaches a critical slope and is slowly driven by random 
input of dropped sand grains. This concept was also applied to the magnetically driven solar corona, which 
releases intermittent amounts of nonpotential magnetic energies during flares (Lu and Hamilton 1991). Many 
applications of SOC processes can be found from geophysics to astrophysics, all the way to financial and 
social systems. The topic of SOC is reviewed in recent reviews, textbooks, and monographs (e.g., Bak 1996; 
Jensen 1998; Turcotte 1999; Charbonneau et al. 2001; Hergarten 2002; Sornette 2004; Aschwanden 2011a; 
Crosby 2011; Pruessner 2012). In this study we go beyond the mere establishment of powerlaw indices, but 
try also to understand the quantitative values of the observed powerlaw slopes and their variation with the 
solar cycle. We interpret the observations in terms of a recently published theoretical model based on a 
statistical fractal-diffusive avalanche model in a slowly-driven self-organized criticality system (Aschwanden 
2012). The application of this model allows us also to provide a quantitative framework for statistical flare 
predictions, which includes also probabilities for the most extreme space weather events. 

The content of this paper is as follows: Section 2 presents the data analysis and results, in Section 3 
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we apply the SOC theoretical model, in Section 4 we compare the results with previous observations, and in 
Section 5 we summarize the conclusions. 

2. DATA ANALYSIS 
2.1. Dataset 

Our basic data input is the 37-year time series of complete years (from 1975 to 2011) of the GOES softer 
X-ray channel in the 1-8 A wavelength range. We choose the softer energy channel (1-8 A = 1.5-12 keV) 
over the harder energy range (0.5-4 A = 3-25 keV) because of its higher sensitivity to smaller flares and 
less data noise. The data are available in the Solar Software (SSW) (Freeland and Handy 1998), originally 
provided by NOAA, and can be read with IDL software (using the procedure RD_GXD.PRO). The data are 
returned with time tick marks U,i — 1,-,N for the 3-s intervals and flux values i = 1,...,N in physical 
units of (Wm~ 2 ), for both wavelength ranges A = 1 — 8 A and A = 0.5 — 4 A. The logarithmic flux values 
are also labeled with letters (A, B, C, M, X-class), which denote the order of magnitude of the peak flux 
on a logarithmic scale (A = 10~ 8 ,B = 10~ 7 ,C = 10" 6 ,M = 10~ 5 ,X = 10~ 4 Wmr' 2 ), subdivided with an 
additional digit (e.g., an X2 class flare has a flux of 2 x 10~ 4 Wm~ 2 ). 

Flare catalogs are issued by NOAA, which provide the peak fluxes, start times, peak times, and end 
times of flare events, but no prcflare background is provided in the flare catalogs, which is important for 
statistical studies of small flares. If no prcflare background is subtracted from the peak flux, the peak flux 
includes not only the flare-associated flux, but also the total soft X-ray flux of all active regions on the 
entire solar disk and beyond the limb, which is larger than the flare-associated flux for A, B and C-class 
flares during active phases of the solar cycle (Wagner 1988; Bornmann 1990; Aschwanden 1994; Veronig et 
al. 2004). Therefore, proper evaluation of the preflare background requires the analysis of the light curves, 
and cannot be obtained from the official NOAA flare catalog alone. Furthermore, the GOES light curves 
contain occasional data gaps (due to Earth occultation, calibration procedures, or drop-outs) and nonsolar 
spikes (of instrumental, terrestrial, or magnetospheric origin) , which need to be removed to avoid false flare 
detections. 



2.2. Automated Flare Detection Algorithm 

After testing various automated flare detection algorithms at various background soft X-ray flux levels, 
which range from the A-level during the solar cycle minimum to the C-level during the solar cycle maximum, 
by optimizing maximum sensitivity of detecting the smallest flares and by minimizing false detections, we 
arrived at the following algorithm: 

1. Rcbinning of data: The intrinsic time resolution of dt = 3 s is rebinncd to time steps of At = 12 s 
(matching the SDO/AIA cadence used in another study) by the median value during each bin with 
length nun = At/dt = 12/3 = 4. Thus, a daily record has riday — 86,400/3 = 28,800 datapoints, 
which is downsampled to a number of nun — 86,400/12 = 7200 bins. 

2. Definition of minimum flare duration: We set a minimum flare duration at a 1-minute time interval, 
i.e., At m in — 60 s, which corresponds to a number of ANbin = 60/12 = 5 bins. 
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3. Definition of noise level and threshold: From the mean and standard deviation of typical GOES fluxes 
during quiescent time periods we find a noise level of f noise ~ 2 x 10 -8 Wm~ 2 , and define a correspond- 
ing threshold level of fthresh = fnoise- The noise level of GOES-8, 9, 10 is constant within less than 10% 



variation ( http://rammb.cira.colostate.edu/l research / calibration_validation_and_visualization / goes Jmage_display / 
noise. asp). 

4. Elimination of data gaps: These are identified by time intervals with a constant "floor" flux value that 
corresponds to the minimum flux of the daily light curve. 

5. Spike removal: Single spikes are detected if the ratio of the maximum flux f max to the minimum flux 
fmin ^ fnoise IS larger than Cfspikes = fmax/fmin > 10 during a time interval of ±Ai m i„ = 60 s. 

6. Smoothing of light curve: We smooth the rebinned light curve with a boxcar of nsm — 21 time bins, 
which corresponds to nsm * 4ANbi n minimum flare durations. 

7. Detection of maxima and minima: We detect now all local maxima and minima of the smoothed light 
curve (consecutively in daily intervals) . The flux maximum times i, represent candidates for flare peak 
times, and flux minimum times and tj+i represent potential flare start and end times. 

8. Detection of flare event: Flare events are defined when they fulfill the following conditions (see Fig. 1): 
(i) The flare starts at a flux minimum time t s = i,-, where a preflare background fsc is defined from the 
median flux in a time interval [t s ~ At m i n , t s ] ; (ii) The flare ends at the first subsequent flux minimum 
time t e = ti + i + k with k > when the flux drops below the level Jbg + fthresh, but before the start 
of a next flare (i.e., no overlapping flare time intervals); (iii) The flare peaks at the highest flux value 
f p at time t p during the intervening time interval t s < t v < t e . (iv) The background-subtracted peak 
flux is F — f p — fsc- Obviously, during solar maximum a lot of flares start on the tail of previous 
larger events, and thus the background flux is higher, which is not the case during solar minimum. The 
GOES non-flare background flux varies similarly to the sunspot number (Wagner 1988). - In addition 
we measure also the steepest flux derivative (df / dt) max during the rise time [t s ,t p ], which can be used 
as a proxi for the hard X-ray peak flux according to the Neupert effect (Dennis and Zarro 1993). 



An example of automated flare detection is shown for 2012-Jan-27 in Fig. 2, when an GOES class XI. 8 
flare occurred. A total of 36 flare events were detected during this day, mostly ranking in the B and C-class 
level. In this example we see that the flux of most of the small B and C-class flares returns to a background 
level of pa 5 x 10~ 7 Wm~ 2 . The duration of larger flares (e.g., C5.5, C.3.3, XI. 8, C3.5) is truncated by the 
start of subsequent flares. Therefore, we expect that the flare duration is generally underestimated at times 
of high flare rates. A more correct flare duration could potentially be derived from a flare decay model that 
extrapolates the beginning of the flare decay phase all the way to the background level. 

We found that this flare detection algorithm is sensitive to the smallest flares recognizable by visual 
inspection, and at the same time to be fairly robust with an estimated false-detection rate of < 10 _2 ...10~ 3 , 
based on visual flare classification. 



2.3. Flare Rate and GOES Duty Cycle 



The monthly rate of detected flares is shown in Fig. 3, for all flares, as well as for the subsets of 
> C1.0, > M1.0, and > A1.0 flares separately. The total number of detected flares larger than a given 
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magnitude decreases approximately by a factor of 10 for each order of magnitude: from 338,661 events of 
all flares to 35,221 > C-class flares, 3986 >M-class flares, to 248 >X-class flares, during a time span of 37 
years. The detected flare rate is also shown on a logarithmic scale (Fig. 3 bottom), which shows better the 
proportionality of detected flares in different magnitudes, even across the minima and maxima of the three 
solar cycles. There are indeed some months during the last extended solar cycle minimum 2008-2009 with 
no detected flares at all, which corroborates also the robustness of our detection algorithm for false events 
by noise coincidence. The number of detected flares per year is listed in Table 1, ranging from a minimum 
of 186 events per year (2008) to a maximum of 18,797 events per year (1979). 

Comparing the number of detected flares with the official GOES flare catalog from NOAA (readable 
with procedure GET.GEV.PRO in SSW/IDL), we find 39,696 flare events reported by NOAA during the 
period of 1991-2011 (the annual number of NOAA events is listed in Table 1, 3rd column). Thus our flare 
detection algorithm detects 9153 events per year in the average, while the NOAA catalog lists 1804 events 
per year, so our algorithm is about 5 times more sensitive. The NOAA flare detection algorithm is defined 
as follows: The event starts when 4 consecutive 1-minute X-ray values have met all three of the following 
conditions: (i) All 4 values are above the Bl threshold; (ii) All 4 values are strictly increasing; (iii) The last 
value is greater than 1.4 times the value that occurred 3 minutes earlier. The peak time is when the flux 
value reaches the next local maximum. The event ends when the current flux reading returns to half of the 



peak value (http://www.ngdc.noaa.gov/stp/solar/solarflares.html). 



We determined the effective GOES duty cycle from the number of datapoints that we obtained from 
reading the GOES database with the SSW standard routine and found an average of 94% ± 4% during the 
years 1978-2011, and a lower fraction of 76% ± 8% during the first 3 years of the GOES series (see Table 
1). The actual duty cycle of GOES may be higher, because our value of the duty cycle is derived from 
the number of readable files, and thus the missing data may also include file reading problems, besides the 
originally missing data (due to data loss, telemetry gaps, calibration procedures, Earth occultation, etc.). 
Anyway, GOES has probably the highest duty cycle from all solar-dedicated space missions, and thus offers 
the most complete record of solar flares over the last three solar cycles. 



2.4. Occurence Frequency Distributions 

The size distribution of many solar flare parameters has been found to be close to a powerlaw function. 
However, often there is a gradual turnover at the lower end of the powerlaw range near the detection threshold. 
In addition, there is sometimes an exponential drop-off apparent at the upper end of the powerlaw range. 
Since the fit of a powerlaw slope can be severely affected by the gradual turnover at the lower end, we fit a 
powerlaw function with an empirical correction term that characterizes the turnover with a transition to a 
constant function at the lower end of the distribution function, 

N{x)dx = NtO- + x/ Xl )~ a dx . (1) 

At the low end the distribution converges to a constant N(x <§; x\) i— > Ni, while the asymptotic limit at the 
upper end turns into a pure powerlaw function, N(x ^> x\) i-> Ni(x/xi)~ a . We call x\ the "lower bound" 
or "turnover value" of the powerlaw range. 

The powerlaw fits of the size distribution functions of the soft X-ray flux are shown for each year from 
1974 to 2012 separately in Fig. 4, yielding powerlaw slopes in the range of ap = 1.7 — 2.3, with a mean and 
standard deviation of ckf — 1-98 ± 0.11. The individual values and uncertainties are also listed for every 
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year in Table 1. The annual fits show in most cases a clear turnover at the lower end, which is well fitted 
by our empirical function (Eq. f ) and renders the value of the powerlaw slope quite stable, independently 
of which bins near the turnover are included or excluded in the fit. The accuracy of the powerlaw slope is 
often hampered in earlier studies due to arbitrary definitions of the fitted range. 

In Fig. 5 we compare the powerlaw fits for the GOES peak fluxes for different phases of the solar cycles. 
We select three different regimes: (i) Years near the solar cycle maximum (defined by periods with more than 
10,000 flare events detected per year, see Table 1, second column, which includes the periods of 1978-1982 
for Cycle 21, 1988-1991 for Cycle 22, 1998-2003 for Cycle 23, and 2011 for Cycle 24), (ii) Years near the 
solar cycle minimum (defined by periods with less than 3000 flare events detected per year, which includes 
the periods of 1975-1976, 1985-1986, 1995-1996, and 2007-2009), and (iii) intermediate time intervals during 
the rise or decay of solar cycles. The powerlaw fits of these three data subsets exhibit a remarkable constant 
powerlaw slope of ccf ~ 2.0 (Fig. 5, top panel), so that the probability distribution functions differ only by 
a variable scaling factor for the occurrence frequency as a function of time. 

In Fig. 5 (middle panel) we show the same statistics for the parameter of the maximum flux-time 
derivative (df / dt) max during the rise time of the flare light curve, which also exhibits the same behavior 
as the soft X-ray peak flux, namely a constant powerlaw slope a d f/ dt ~ 2.0, with a variable scaling factor 
for the occurrence frequency as a function of time. The powerlaw slopes a d f / dt are also tabulated for each 
year in Table 1. If we would apply the Neupert effect, the time derivative (df /dt) of the soft X-ray flux 
should be a good proxi of the (smoothed) hard X-ray peak flux, which predicts then a powerlaw slope of 
a P p ed — a d f/ d t = 2.0. In reality, however, hard X-ray peak rates were found to have a powerlaw slope of 
ap = 1.73 ± 0.07 (Aschwanden 2011b). This discrepancy results from a combination of two effects: (i) 
The simplest formulation of the Neupert effect in terms of an integral (Eq. 21) is oversimplified as we will 
see below (Section 3.2), and (ii) the time derivative measured from the smoothed soft X-ray time profile 
represents an underestimate of the true time derivative of the unsmoothed light curve. 

Finally, we plot also the frequency distribution of flare rise times t r i se , which is probably a good proxy 
for the flare duration T of effective energy release, as measured from non-thermal X-rays light curves and 
expected from the Neupert effect. We plot the size distributions for the same three epochs of the solar cycle, 
but find a completely different behavior (Fig. 5, bottom panel). The powerlaw slopes are found to be steepest 
at the solar cycle maximum with an average slope of ap ~ 3.2, while they become flatter with a value of 
«t ~ 2.3 at the solar minimum. The values for different years vary from ay = 1-75 (for the year 2008 with 
the most extreme solar minimum), up to olt — 5.18 (for the year 2002, which represented the peak of the 
last solar Cycle 23), as it can be seen from the list in Table 1. 

2.5. Solar Cycle Variatiability 

Apparently, some soft X-ray parameters vary during the solar cycle, while others do not. For an overview 
of this time-dependent behavior we plotted the various parameters as a function of time (over the last 37 
years) in Fig. 6. The most invariant parameters are the powerlaw slope ap ~ 2.0 of the soft X-ray peak flux 
and the slope ct d f/ d f « 2.0 of the time derivative, having almost identical values. The largest deviation from 
a constant value is seen for the year 2009, the year with the second- lowest flare rate (Nfi are — 528), and 
thus represents a statistical fluctuation due to small- number statistics. 

The time variation of the (detected) soft X-ray flare rate is shown in the bottom panel of Fig. 6 (hatched 
with a grey area). Comparing with this solar cycle variation of the flare rate, we see clearly that the other 
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parameters plotted in Fig. 6 all exhibit a correlated time variation, which applies to the powerlaw slope olt 
of the flare rise time T (Fig. 5, 3rd panel), the turnover value Pi of the peak count rate P (corresponding 
to the symbol x\ in Eq. 1) (Fig. 5, 4th panel), the turnover value (df /dt)\ of the time derivative (df/dt) 
(Fig. 5, 5th panel), and the turnover value T\ of the rise time T (Fig. 5, 6th panel). In all these parameters, 
the turnover value seems to fluctuate in synchronization with the flare rate. This implies that the degree 
of undersampling of weak flares becomes more severe for times with higher flare rates, especially during 
the solar cycle maxima. This behavior is apparently a feature of our automated flare detection algorithm, 
which detects a smaller number of flares during periods of high activity, either because the flux contrast for 
small flares is relatively weaker on top of a light curve of a large flare, compared with the background noise 
level during quiescent times. In addition, the rule that each flare is separated in time and does not allow 
for temporal overlap (of rise times), may suppress the recording of smaller flares during the evolution of a 
larger flare. Actually, we will demonstrate in the following theoretical section that this effect can be simply 
modeled and predicted as a function of the instantaneous flare rate (Section 3.3). 

3. THEORETICAL MODEL 

We apply now the statistical fractal-diffusive (FD-SOC) avalanche model of a slowly-driven self-organized 
criticality (SOC) system, which is derived in general form for all 3 Euclidean dimensions S = 1,2,3 and 
tested with cellular automaton simulations in a recent paper (Aschwanden 2012). For our application to soft 
X-ray data from solar flares here we consider only the 3D case (S — 3) and add a simple physical model of 
the soft X-ray (thermal) radiation mechanism (Section 3.2). 

3.1. The Fractal-Diffusive Avalanche SOC Model 

The fractal-diffusive SOC model is a universal (physics-free) analytical model that describes the statis- 
tical time evolution and occurrence frequency distribution function of SOC processes. It is based on four 
fundamental assumptions: (1) A SOC avalanche grows spatially like a diffusive process; (2) The spatial 
volume of the instantaneous energy dissipation rate is fractal; (3) The time-averaged fractal dimension is 
the mean of the minimum dimension Ds.mhi ~ 1 (for a sparse SOC avalanche) and the maximum dimension 
Ds,max = S (given by the Euclidean space); and (4) The occurrence frequency distribution of length scales 
is reciprocal to the size L of spatial scales, i.e., N(L) oc L~ s in Euclidean space with dimension S. We will 
discuss these assumptions in more detail in the following. 

The first assumption of a diffusive process is based on numerical simulations of cellular automaton 
models. A SOC avalanche propagates in a cellular automaton model by next-neighbor interactions in a 
critical state, where energy dissipation propagates only to the next neighbor cells (in a S-dimensional lattice 
grid) that are above a critical threshold. This mathematical rule that describes the entire dynamics and 
evolution of a SOC avalanche is very simple for a single time step, but leads to extremely complex spatial 
patterns after a finite number of time steps. For a visualization of a large number of such complex spatial 
patterns generated by a simple iterative mathematical redistribution rule see, for instance, the book "A 
New Kind of Science" by Wolfram (2002). The complexity of these spatial patterns can fortunately be 
characterized with a single number, the fractal dimension D$. If one monitors the time evolution of a spatial 
pattern of a SOC avalanche in a cellular automaton model, one finds that the length scale x(t) evolves with 
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time approximately with a diffusive scaling (Section 2.1 in Aschwanden 2012), 

x(t) oc t 1 ' 2 , (2) 

which leads to a scaling law between the avalanche sizes L = x(t = T) and time durations T of SOC 
avalanches, 

L oc T 1 ' 2 . (3) 

The second assumption of a fractal pattern of the instantaneous energy rate is also based on tests with 
cellular automaton simulations (see Fig. 2 in Aschwanden 2012). The fractal dimension is essentially a 
simplified parameter that describes the "micro-roughness" , "graininess" , or inhomogencity of critical nodes 
in a lattice grid in the state of self-organized criticality. Of course, such a single number is a gross over- 
simplification of a complex system with a large number of degrees of freedom, but the numerical simulations 
confirm that it reproduces the correct scaling law between the instantaneous energy dissipation volume Vs it) 
and the spatial scale x(t) of SOC avalanches, 

V s (t)^x Ds , (4) 

which leads also to a statistical scaling law between avalanche volumes V and spatial scales L or durations 
T (with Eq. 3) of SOC avalanches, 

V s oc L Ds oc T Ds ' 2 . (5) 



The third assumption of the mean fractal dimension has also been confirmed by numerical simulations 
of cellular automaton SOC processes in all three dimensions S = 1,2,3 (Aschwanden 2012), but it can 
also be understood by the following plausibility argument. The sparsest SOC avalanche that propagates by 
next-neighbor interactions is the one the spreads only in one spatial dimension, and thus yields an estimate 
of the minimum fractal dimension of Ds, m in ~ 1, while the largest SOC avalanche is almost space- filling and 
has a volume that scales with the Euclidean dimension, Ds, m ax — S. Combining these two extreme cases, 
we can estimate a time-averaged fractal dimension from the arithmetic mean, 

+ D S .max 

< D S — ^— = —, (6) 

which yields a mean fractal dimension of < D$ > = (1 + 3)/2 = 2.0 for the 3D case (S = 3). 

The last assumption of the size distribution is a probability argument. The system size L sys of a 
SOC system represents an upper limit of spatial scales L for SOC avalanches, i.e., L < L sys . For the 
3D-case, the volumes V of individual avalanches are also bound by the volume V sys of the system size, i.e., 
V = L 3 < V sys = Ll ys . If the entire system is in a critical state, SOC avalanches can be produced everywhere 
in the system, and the probability N(L) for a fixed avalanche size L with volume V is simply reciprocal to 
the size, i.e., 

N(V) oc oc V^ 1 , (7) 



which is equivalent to 



L 3 

N(L) oc -p— oc L~ 3 . (8) 



Based on this model we simulate now an example of a time evolution of SOC parameters for the 3D 
case, as shown in Fig. 7. We start with a mean fractal dimension < D 3 >= 2.0, and simulate the time 
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evolution of the fractal dimension D^{t) by simulating fluctuations with a relatively small amplitude (with 
a standard deviation of <td = O.f 5) by a random generator (Fig. 7 top panel). The time evolution of the 
instantaneous energy dissipation rate f(t) in classical SOC systems (with a constant mean dissipated energy 
quantum < AE > per unstable lattice node) is then proportional to the instantaneous volume Vs(t), which 
yields the following function (with Eq. 4), 

/(*) = oc< AE > V s (t) =< AE > x(t) Ds =< AE > t Ds ' 2 . (9) 

In the 3D case with D$ — 2, we expect than the proportionality f(t) oc t 1-0 , as shown in Fig. 7 (second 
panel). 

The statistical peak value p(t) of the energy dissipation rate after time t can be estimated from the largest 
possible avalanches, which have an almost space-filling dimension D$ < S, and thus would be expected to 
scale as p(t) oc t s / 2 oc t 1 - 5 (indicated with a dashed curve in Fig. 7, second panel. 

The evolution of the total dissipated energy e(t) after time t is simply the time integral, for which we 
expect 

e( t ) = f ^lldr oc f t Ds ' 2 oc t< 1+Ds M , (10) 
Jo dT Jo 

which yields the function e(t) oc t 2 for the 3D case (Fig. 7, third panel). These time evolutions apply to 
every classical SOC model and can be derived from purely statistical probability arguments, and thus are 
physics-free. For application to real- world data we have to understand the physical nature of the observables, 
before we can relate the measured distributions to the statistical SOC theory, which we undertake in the 
next Section. 



3.2. SOC Application to Solar Flare Soft X-ray Data 

SOC theory has been applied to different wavelength regimes of solar flare observations, such as to 
gamma-rays, hard X-rays, soft X-rays, and extreme ultra-violet (EUV). A comprehensive review of such 
studies is given in Section 7 of Aschwandcn (2011a). However, since each wavelength range represents 
a different physical radiation mechanism, we have to combine now the physics of observables with the 
(physics-free) SOC statistics. 

Soft X-ray emission during solar flares is generally believed to result from thermal free-free and free- 
bound radiation of plasma that is heated in the chromosphere by precipitation of non-thermal electrons 
and ions, and which subsequently flows up into coronal flare (or post-flare) loops, a process called "the 
chromosphcric evaporation process" (for a review see, e.g., Aschwanden 2004). Therefore, we can consider the 
flare-driven chromospheric heating rate as the instantaneous energy dissipation process of a SOC avalanche, 
as shown in the simulated function f(t) in Fig. 7 (second panel). The heated plasma, while it fills the 
coronal flare loops, loses energy by thermal conduction and by radiation of soft X-ray and EUV photons, 
which generally can be characterized by an exponential decay function after an impulsive heating spike. In 
Fig. 7 (bottom) we mimic such a soft X-ray radiation light curve by convolving the instantaneous energy 
dissipation rate f(t) (Fig. 7, second panel) with an exponentially decaying radiation function (with an 
c- folding time constant of Td ecay ), 



fsxr(t)= / f(t)exp[-(t-t')/Tdecay]dt' , 
J — oo 



(11) 
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which shows also a time dependence that follows approximately 

/a*r(t) OC /(f) CX t 1 " , (12) 

because the convolution with an exponential function with a constant e-folding time constant acts like a 
constant multiplier. In the limit of infinitely long decay times (rdecay °°), our convolution function 
(Eq. 11) turns into a time integral of the heating function f(t), which is also known as Neupcrt effect 
(Dennis and Zarro 1993; Dennis et al. 2003), where the heating function is identified with the non-thermal 
hard X-ray emission and the time integral with the soft X-ray emission. However, such an approximation 
would predict different values for the powerlaw slopes of hard X-ray and soft X-ray peak fluxes, which was 
found to contradict the data (Lee et al. 1995). Hence, we consider the formulation of the Neupcrt effect in 
terms of a convolution with a finite cooling time (Eq. 11) as a more accurate representation than the time 
intcrgral formulation (with an infinite cooling time). 

We can now calculate the occurrence frequency distributions. The total duration T of energy release 
corresponds essentially to the rise time t r i se of the soft X-ray flux, because the decay phase of a soft X-ray 
flare light curve is dominated by conductive and radiative loss, rather than by continued heating input. Thus, 
starting from the size distribution of length scales, N(L) oc L~ s (Eq. 8), we can derive the size distribution 
of rise times T = t r i se by substituting the variable T for L in the distribution N(L), using the diffusive 
scaling L(T) oc T 1 ' 2 (Eq. 3), with the derivative dL/dT = T" 1 / 2 , 

dT cx T- [{1+s),2] dT = T- 2 dT , (13) 

yielding N(T) oc T~ 2 for 5 = 3. Subsequently we can derive the size distribution N(F) of the energy 
dissipation rate F = f(t = T), using the relationship F(T) oc T Ds / 2 (Eq. 9), 

dF oc F -H+(s-i)/d s ] = F -2 dp (u) 

Thanks to the proportionality between the functions /(f) and f sxr (t) (Eq. 12), we have also a proportionality 
between the statistical expectation values F = f(t = f) and F sxr — f sxr (t), and thus a size distribution with 
an identical powerlaw index, 

N \F sxr ^dF sxr oc F s J- r ^ )/ ] F sxr dF sxr . (1^) 

Thus, our theory predicts a powerlaw slope of a_F = 2.0 for the soft X-ray fluxes (Eq. 15), and an identical 
slope of olt = 2.0 for the rise time t r ise — T (Eq. 13). As we can see in the observations shown in Fig. 6, 
the powerlaw slopes of the soft X-ray fluxes has an annually averaged powerlaw slope of qj? = 1.98 ± 0.11, 
which is moreover invariant over a time span of 37 years, and thus is fully consistent with our theoretical 
prediction of Of = 2.0. The powerlaw slopes of soft X-ray rise times have a minimum value of ax > 2.0 
during the solar cycle minima (Fig. 6, third panel), but deviate up to a T < 5.0 during solar cycle maxima. 
Apparently there is a solar cycle effect that modifies the theoretically predicted value, which we interpret as 
a measurement bias in the next Section. 



3.3. Solar Cycle Dependence of Flare Pile-Up Bias 

Our theory predicts an occurrence frequency distribution of N(T) oc T~ 2 for flare durations (or rise 
times in soft X-ray light curves) in the slowly-driven limit where subsequent avalanches do not overlap in 



N(T)dT = N(L[T}) 



dL 
dT 



N(F)dF = N(T[F]) g 
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time. Also, our automated flare detection algorithm does not allow for overlapping flare time intervals. 
However, it is conceivable that the flare rate is so high during the most active times of the solar cycle that 
multiple flares overlap each other, which violates the slowly-driven condition and leads to an underestimate 
of the flare durations (because of our rule that a flare has to end before the next flare starts), especially 
during large and long flares, a bias that leads to a steepening of the powcrlaw slope. Indeed we observe 
powerlaws that are steeper up to a value of ojt = 5.18 ± 0.55 (in the year 2002 near the maximum of the 
last solar cycle). Thus we try to estimate this effect, which we call "flare pile- up bias". 

The average waiting time between two subsequent flares per year (with a duration of 365.25 x 86, 400 s 
= 3.2 x 10 7 s) is 

3.2xl0 7 s 

< At wait >= — , (16) 

JV flare 

which varies between < At wait > = 3.2 x 10 7 s/18, 797 w 1700 s « 28 min (in the year 1979) to < At wait >= 
3.2 x 10 7 s/186 w 170,000 s « 2 days (in the year 2008). The latter year is sufficiently low to fulfill the 
slowly-driven condition, because the flares are so sparse that they never overlap. However, the former year 
with an average waiting time of 1.7 x 10 3 s is significantly shorter than the longest measured flare rise time 
(trise ~ 10 4 ' 5 s) and thus underestimates the longer flare durations. We estimate the effect of this flare pile-up 
bias on the measurement of the powerlaw slope, which is generally defined as, 

a T = -»1 (17) 

log (T 2 /ro ' [U) 

with T\ and T 2 the lower and upper bound of the powerlaw range, and Ni and N 2 the corresponding number 
of flares at these boundaries. In the case of flare pile-ups the upper limit T 2 is limited by a time scale that 
approximately corresponds to the average waiting time < At wa it >, modifying the observed powerlaw slope 
a^ bs to 

a o bs = logjNM) 
T log (< At wait > /Ti) ' V ; 

for < At wait >< T 2 . From years with sparse flare rates fulfilling the slowly-driven condition we find 
log(N 2 /N 1 ) « 5.2, log (Ta/Ti) w 2.6, yielding a$ s w 2.0. Combining Eq. (16) and (18) we obtain a time- 
dependent powerlaw index that depends on the annual flare rate Nfi are (t) as, 

1 j ~ log(3.2xl0 7 /7V /iare WT 1 ) ' (i9j 

if Nfi are > 3.2 x 10 7 s/T 2 - We calculate this time-dependent powerlaw index s (t) for all years (1975-2011) 
and display the expected curve in Fig. 6 (3rd panel), which closely mimics the powerlaw rate variation a^ s {t) 
inferred from the observations. Thus, we conclude that the correct powerlaw slope of cut = 2.0 predicted by 
SOC theory can be measured during solar cycle minima, while flare pile-up causes a steeper powerlaw slope 
during solar cycle maxima. 



4. DISCUSSION 
4.1. Comparison with Previous Statistics 



The size distributions of soft X-ray peak fluxes of solar flares have been reported from various spacecraft 
data, with powerlaw indices in the range of ap w 1.6 — 2.1 (Hudson et al. 1969, Drake et al. 1971, Shimizu 
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1995, Lee et al. 1995, Feldman et al. 1997, Shimojo and Shibata 1999, Veronig et al. 2002a,b, Yashiro et al. 
2006). A compilation of these studies is listed in Table 2. Let us quickly review these results, particularly 
under the aspect of compatibility with our new results. 

Probably the first size distribution of soft X-ray fluxes from solar flares was reported for 177 events 
observed with OSO-3 by Hudson et al. (1969). The cumulative size distribution shows a powerlaw slope of 
ftp w 0.85 and a gradual turnover at the lower end, which corresponds to ap — f3p + 1 m 1.85 and is roughly 
consistent with SOC theory (ap = 2), given the small- number statistics. 

A larger statistics of 3140 flare events was gathered with the Explorer spacecraft (Drake et al. 1971), 
finding a powerlaw slope of ap — 1.84 ±0.02 if the fit extended to the third-lowest bin, and a flatter value of 
ap = 1.66 ± 0.02 if the fit extended to the second-lowest bin. Clearly, the turnover at low values flattens the 
powerlaw slope as expected, a systematic error that is not included in the error bars. The same study finds 
also a powerlaw slope of ap — 1-44 ± 0.01 for the fluence, and thus is largely consistent with SOC theory 
(a P = 2.0, a E = 1.5). 

Shimizu (1995) sampled active region brightenings with Yohkoh, which we consider as "miniature flares" . 
Measuring the peak soft X-ray intensities averaged over 8x8 macropixels, he found a powerlaw slope 
of ap = 1.80, and ap = 1.68 for 16 x 16 macropixels, while the powerlaw of energies was found to be 
aE ~ 1.5 — 1.6, which is largely consistent with SOC theory (ap = 2.0, = 1.5). Another study on Yohkoh 
active region brigthenings reported a less accurate value of ap = 1.7 ± 0.4 (Shimojo and Shibata 1999), due 
to the relatively small sample of 92 events. 

Lee et al. (1995) reported powerlaw slopes of ap — 1.86 for soft X-ray peak fluxes observed with 
BCS/SMM, and aslope of a P = 1.86 for a sample of 4356 GOES events. A similar value of a P = 1.88 ±0.21 
was reported by Feldman et al. (1997), all within about 5% of the theoretically predicted value. 

A larger dataset of GOES flares during the years of 1986-2000 was analyzed by Veronig et al. (2002a), 
yielding a powerlaw index of ap = 1.98 ± 0.08 after preflare background subtraction, which is perfectly 
consistent with SOC theory and the values obtained in this study (ap = 1.98 ± 0.11). In a related study 
with 49,409 GOES flares, where no preflare background was subtracted, a value of ap = 2. 11 ±0.13 was found 
(Veronig et al. 2002b). Since the powerlaw range covered only medium-size to large flares, the background 
subtraction that is important for small events did not spoil the powerlaw fit. The powerlaw fit of soft X-ray 
flare durations was found to be ap = 2.93 ± 0.12. This is not inconsistent with our findings, because this 
powerlaw index was found to vary in the range of ax ~ 2 — 5 from the minimum to the maximum of the 
solar cycle, which we interpret as flare pile-up bias. 

Very similar values were reported for GOES flares during the period of 1996-2005 (Yashiro et al. 2006), 
with ap = 2.16 ± 0.03 for all flares, and ap = 1.98 ± 0.05 for flares with CMEs. The complementary subset 
of flares without CMEs have a higher value of ap = 2.52 ± 0.03, which can be understood as a selection 
effect that preferentially excludes the larger flares, and thus is expected to produce a steeper powerlaw slope. 
The flare durations showed powerlaw slopes in the range of ap ~ 2.5 — 3.2 for the same subsets, which can 
also be understood in terms of the flare pile-up bias. 

In summary, all previous measurements are consistent with the theoretically predicted powerlaw slope 
of ap = 2.0, if we take the uncertainties of small-number statistics and the turnover at the lower bound of 
the powerlaw range into account. The distribution of flare durations was generally reported to be higher 
(ap ~ 2.5 — 3.2) than the theoretically predicted value ap = 2.0, but can be satisfactorily understood in 
terms of the flare pile-up bias effect. 
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4.2. 



Coronal Heating by Nanoflares 



If a powerlaw distribution of flare energies E is derived, the critical slope of aE,crit — 2.0 decides whether 
the integral of the differential flare energy distribution diverges at the lower or upper end of the powerlaw 
range (Hudson et al. 1991). Our SOC model, which is increasingly supported by previous observations 
(Aschwanden 2011b, c) as well as with the new analysis presented here, predicts a powerlaw slope of ole = 1-5 
for total time-integrated energies, which is also consistent with observed occurrence frequency distributions 
in hard X-rays (e.g., ue = 1.53 ±0.02; Crosby et al. 1993). Based on this agreement between the theoretical 
and observational results there is neither a theoretical prediction nor observational evidence for a powerlaw 
slope ole steeper than the critical value of aE,crit = 2. In addition, the theoretical SOC model predicts a 
powerlaw slope of ap = 2.0 for peak energy dissipation rates, which requires an energy slope of as = 1.5 
to be consistent with SOC theory. This has also been confirmed in this study, which further corroborates 
the theoretical SOC model. Based on this argument we conclude that nanoflares have only a negligible 
contribution to coronal heating. 



What do we learn about the physics of solar flares by applying a statistical SOC model? The classical 
cellular automaton model mimics the complex spatial and temporal patterns that result in a nonlinear 
energy dissipation process, without involving a physical mechanism. However, it is thought that a solar flare 
consists of a fragmented energy release process (Benz 1985), which is manifested in a highly intermittent 
time structure of dccimctric radio emission or nonthermal hard X-ray emission (with time structures down 
to milliseconds). Numerical MHD simulations can reproduce such a intermittent time evolution in chain 
reactions of tearing-mode instabilities and magnetic island formation (e.g., Sturrock 1966; LeBoef et al. 
1982; Tajima et al. 1987; Kliem 1990, 1995; Karpen et al. 1995; Drake et al. 2006a,b), leading to fractal 
current sheets (Shibata and Tanuma 2001). What the SOC model tells us are the scaling laws between the 
energy dissipation rate, total dissipated energy, spatial scales, and flare durations, since the power indexes of 
the correlated parameters are linked to the powerlaw distributions of these parameters. While SOC theory 
provides the statistical framework of basic spatio-temporal parameters (such as time scales, length scales, 
and fractal volume scaling), the physics comes into the problem by connecting the observables (i.e., flux 
intensities at different wavelengths) to the geometric SOC parameters, which involves physical modeling in 
terms of electron densities, electron temperatures, thermal, and nonthermal energies. 

In this study we found that the powerlaw slopes are consistent with the assumption that the emitted soft 
X-ray flux during a solar flare is proportional to the fractal flare volume (Eq. 9), i.e., f sxr {t) oc f(t) oc Vs(t). 
Given the fact that soft X-ray emission consists of thermal free-free emission, which in the optically-thin 
limit has a proportionality of the emission measure to the squared density in a volume element cell dV, 



the emissivity per volume cell is expected to depend on the electron density. Physical modeling is required 
to quantify the soft X-ray emissivity and emission measure per volume element in a given wavelength range 
as a function of the volumetric heating rate, conductive, and radiative loss rate (e.g., Metcalf and Fisher 
1995; Phillips and Feldman 1995; Garcia 1998, 2001; Aschwanden and Alexander 2001; Battaglia et al. 2005; 
Sylwester et al. 1995; Veronig and Brown 2004; Veronig et al. 2005; Yamamoto and Sakurai 2010). In sum- 
mary, SOC theory represents an important statistical tool to link physical modeling to observed correlations 



4.3. 



Physical Aspects of the SOC Model 




(20) 
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(scaling laws) and size distributions. 



4.4. 



Soft X-ray versus Hard X-ray Flare Statistics 



A compilation of hard X-ray flare statistics over the last three solar cycles has been conducted in a 
recent study (Aschwanden 2011b and references therein). The following powerlaw slopes have been found for 
the hard X-ray parameters: ap — 1.73 ± 0.07 for the hard X-ray peak flux P, (Xe — 1-62 ± 0.12 for the total 
time-integrated hard X-ray flux (or fiuence) E, ar — 1-99 ± 0.35 for the hard X-ray flare durations T . These 
results agree remarkably well with ouf FD-SOC model, which predicts ap = 5/3 « 1.67, ole = 3/2 = 1.50, 
and olt = 2.00 for the 3D case (Aschwanden 2012; Table 1 therein). Comparing hard X-ray with soft X-ray 
flare statistics, we have to be aware that the peak fluxes in these two wavelengths are not equivalent. The 
peak of the hard X-ray flux light curve fhxr{t) represents a maximum fluctuation of the instantaneous energy 
dissipation rate (if the temporal fluctuations are fully resolved in time), while the peak flux in soft X-rays 
represents the maximum value of the smoothed energy dissipation rate, because the intermittent and spiky 
instantaneous energy dissipation rate is convolved with a cooling time, which has a smoothing effect on the 
light curve, as demonstrated in Fig. 7. Therefore, our FD-SOC model predicts a steeper slope cxf = 2.0 for 
the soft X-ray peak flux distribution than the slope ap = 1.67 for the hard X-ray peak flux distribution, 
which is indeed confirmed with the observed data. 

In an earlier study (Lee et al. 1995) it was suspected that the size distribution of the soft peak flux 
(with powerlaw slope oif) should correspond to the size distribution of the (time- integrated) hard X-ray 
fiuence ae, according to the Neupert effect, which predicts 



but the data were found to contract this expectation (Lee et al. 1995). This failure corroborates that the 
integral formulation (Eq. 21) of the Neupert effect is an oversimplified approximation and does not hold on 
a statistical basis. A more accurate approximation is the formulation of the Neupert effect in terms of a 
convolution of the hard X-ray flux with a finite e-folding cooling time, as expressed with Eq. (11). Such 
a relationship is consistent with the predictions of our FD-SOC model, which predicts a powerlaw slope of 
ap = 2.0 for the soft X-ray peak flux and a powerlaw slope of = 1.5 for the hard X-ray fiuence (while the 
Neupert effect would predict both to have an identical value) . The relationship between soft X-ray and hard 
X-ray fluxes has also been quantified in some statistical studies (e.g., Battaglia and Benz 2006; McTiernan 
2009; Falewicz et al. 2009). 

The powerlaw slope ap of the hard X-ray peak flux was found to vary sligthly with the solar cycle, 
between a minimum value of ap = 1.62 ± 0.02 in the decay phase of a solar cycle to ap = 1.79 ± 0.02 during 
the minimum or rise of the solar cycle (Aschwanden 2012a). However, this variation is about of the same 
order as the differences between the three different instruments (SMM, CGRO, RHESSI) and three different 
flare detection methods used for the statistics. Since we find that some soft X-ray parameters are invariant 
during the solar cycle, while other parameters vary due to the flare pile-up bias, it is not clear at this point 
whether the variation of hard X-ray powerlaws is also caused by the flare pile-up bias, or if it is related to 
intrinsic variations of the self-organized criticality conditions (Aschwanden 2011c). 




(21) 
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4.5. Solar Flare Predictability 

The occurrence frequency distributions we obtained from the GOES soft X-ray fluxes are also called 
probability distribution functions (PDF) and can be used for statistical flare prediction. The fact that we 
established here the invariance of the powerlaw slope, makes the statistical prediction very easy, because 
the number of flares observed or to be predicted depends only on one single time-dependent parameter, the 
scaling factor Ni, 

N(P)dP = m(t) (j^j P dP = N 1 (t)^ dP, (22) 

where Pi is an absolute constant, say Pi — 10~ 6 Wm" 2 for the C-class level, and iVi(i) is the number of 
flares observed at this level in the interval dP. Forecasting the number of flares that are larger than a given 
threshold level, are given by the cumulative PDF, 

N cum (> P) = £ Ni(t) (-0 dP = Ni(t) (J^j . (23) 

The time- varying function Ni(t) can be measured at any level Pi within the powerlaw range. This time- 
varying function iVi (t) can be broken down for different active regions and be turned on and off individually 
for each active region in a prediction model, depending on their appearance on the east-side of the Sun or 
disappearance on the west-side. Thus, the accuracy of forecasting depends only on the temporal extrapolation 
of the flaring rate Ni(t) per active region, which can also be trended from statistical variations. 

Predicting the largest event in history or future is a delicate issue. In principle we can take the mean 
of the largest events during each solar cycle as a reasonable guess for the next solar cycle, which yields 
something like an X20 GOES class. The solar dynamo resets the magnetic field after every 11-year cycle, so 
that no nonpotential magnetic field can be stored on longer terms, potentially leading to larger giant flares. 
What the SOC theory tells us in addition is the scaling of the flare energy with the length scale, which simply 
scales with the fractal volume, E oc Vs oc L Ds . Once the size scale L is measured for the largest flares, we 
can scale the maximum possible flare energy by the maximum scale size L max , which probably corresponds 
to the maximum active region size Lar- Active regions grow and wane in size, LAR,{t), which can be used 
to estimate the maximum flare rate Ni{t) per active region, and this way be used in the prediction of the 
cumulative PDF of flares above some size (Eq. 23). 



5. CONCLUSIONS 

In order to obtain a deeper physical understanding of nonlinear dissipative systems governed by self- 
organized criticality (SOC) we need large statistics of SOC parameters to accurately quantify their statistical 
probability distribution functions (PSF) and the scaling laws between SOC parameters. For this purpose 
we analyzed the largest available, uniformly sampled dataset of solar flare light curves from the GOES 
spacecraft series over the last 37 years (1974-2012). Analyzing this dataset and applying the fractal-diffusive 
self-organized criticality model (FD-SOC) we arrived at the following major conclusions: 

1. With an automated flare detection algorithm applied to the GOES 1-8 A light curves we detect a total 
of 338,661 flare events during the epoch of 1975-2011 that includes all full years with GOES records. 
Our algorithm detects about 5 times more solar flares (with a noise level at the GOES A2-class level) 
than the official GOES flare catalog issued by NOAA, which covers the epoch of 1991-2012. The duty 
cycle of GOES is found to be 94%±4% during the years 1978-2011. 
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2. The occurrence frequency distributions of GOES flare parameters can be characterized by a powerlaw 
function with a gradual turnover at the lower end, which can be characterized and fitted with the 
function N(x) = Ni(l + xjx\)~ a , in order to obtain an accurate powerlaw slope corrected for the 
turnover. We find the following mean powerlaw slopes: ap = 1.98 ± 0.11 for the soft X-ray peak flux 
F, ap = 2.97 ± 0.71 for the soft X-ray rise time T, and etdf/dt — 2-01 ± 0.12 for the steepest time 
derivative during the soft X-ray rise time. 

3. The powerlaw slopes of the (background-subtracted) peak flux F and of the time derivative (df/dt) 
are found to be invariant during the last three solar cycles with a variation of less than 5%. The 
powerlaw slope of the soft X-ray rise time T and the turnover values (F\,df /dt\,T\) are all found to 
vary systematically with the solar cycle and can be modeled with the flare pile- up bias effect, which 
causes a loss of small flares and an underestimate of the flare duration (defined by the rise time here) 
during times of high flare rates. The flare rate of years during the solar minimum, however, are not 
affected by the flare pile-up bias, during which we measure a powerlaw slope of «t = 2.02 ± 0.04 for 
the soft X-ray rise time T. 

4. The fractal-diffusive self-organized criticality (FD-SOC) model (Aschwanden 2011a) predicts a fractal 
dimension of Z) 3 = 2.0 for flare volumes and powerlaw distributions with slopes of ap = 2.0 for the 
(smoothed) peak flux F and a T = 2.0 for the flare duration, which are both consistent with the 
observations (taking the flare pile- up effect into account). Most previous studies are also consistent 
with our results and the FD-SOC model (if we take flare pile- up effects into account), except those 
with small-number statistics or no background subtraction. 

5. The consistency of the observations with the theoretical (FD-SOC) model implies also a powerlaw 
distribution for the time-integrated dissipated flare energy with a powerlaw slope of <xe = 1-5, which 
rules out that nanoflares contribute significantly to coronal heating (which would require a powerlaw 
slope of up > ct crit = 2.0.) 

6. The theoretical FD-SOC model predicts the powerlaw distributions of the instantaneous energy dis- 
sipation rate F, the peak dissipation rate P, the time-integrated energy release volume E, the time 
duration T, the fractal dimension D$ of the avalanches volume Vs, and scaling laws between these 
spatial and temporal parameters of SOC avalanches. The prediction of the size distributions of ob- 
servables, such as the soft X-ray flux F sxr or hard X-ray flux Fhxr requires a physical model of the 
emissivity per volume element, which entails thermal and non-thermal emission. For soft X-ray data 
we find that the soft X-ray intensity Fsxr is proportional to the instantaneous energy dissipation rate 
F, which implies a balance between the plasma heating rate and the radiative loss rate. 

7. The hard X-ray peak fluxes have a powerlaw slope ap « 1.7 that is different from the powerlaw slope 
ap = 2.0 of the soft X-ray peak flux, but both values are consistent with the FD-SOC theory. The 
difference can be understood by the fact that the hard X-ray peak flux represents the most extreme 
fluctuation of the energy dissipation rate during a flare, while the soft X-ray peak flux represents 
a convolution function with an e-folding cooling time, which smoothes out the fluctuations of the 
intermittent energy dissipation rate. The integral formulation of the Neupcrt effect appears to be an 
oversimplification that is not consistent with the data, but a formulation in terms of a convolution of 
the hard X-ray flux with an e-folding cooling time seems to be more realistic and is consistent with 
the FD-SOC model. 

8. The invariance of the soft X-ray peak flux powerlaw slope (ap = 2) during the solar cycles simplifies 
the statistical prediction of solar flares, since the statistical expectation value is then proportional for 
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every flux level. The cumulative occurence frequency distribution has a slope of f3p = ctp — 1 = 1, 
which states that the probability of a flare greater than a flux level F drops reciprocally with the flux 
level F. An X-class flare is 100 times less likely than a C-class flare. The largest flares in history or 
future are unlikely to surpass the largest flares observed hitherto, based on the fact that the FD-SOC 
model predicts a dependence of the released energy on the flare size, which is limited by the largest 
system size, i.e., the largest active region in the case of solar flares. 

If the reader goes back to the eight questions at the beginning of the Introduction, he/she will find that 
we arrived at a quantitative answer to almost all questions. Nevertheless, there are more open questions: 
How can we understand the solar cycle variation of the hard X-ray powerlaw slope? What is the physical 
scaling between the (fractal) avalanche volume and the emissivity in different wavelengths (soft X-rays, hard 
X-rays, gamma-rays, EUV)? How do the scaling laws between different forms of flare energies relate to each 
other (thermal, non-thermal, kinetic, magnetic energy)? Does a more extended flare energy distribution 
spanning over 8 orders of magnitude from the largest flare to the smallest detectable nanoflares exhibit a 
single powerlaw slope? Can we statistically verify the fractal spatial geometry and scaling laws with imaging 
observations of solar flares? Future studies with multi- wavelength data from AIA/SDO are expected to shed 
more light on these problems. 

The first author acknowledges helpful discussions with Karel Schrijver and Greg Slater. We thank 
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Table 1. Powcrlaw slopes of occurrence frequency distributions detected with GOES per year. 
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Table 1 — Continued 



Year 
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Powerlaw slope 
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Table 2. Frequency distributions measured from solar flares in soft X-rays. References: 1 ) Hudson et al. 
(1969); 2 ) Drake et al. (1971); 3 ) Shimizu (1995); 4 ) Lee ct al. (1995); 5 ) Feldman et al. (1997); 6 ) Shimojo 

and Shibata (1999); 7 ) Vcronig et al. (2002b); 8 ) Veronig et al. (2002a); 9 ) Yashiro et al. (2006); 10 ) This 
work; n ) Aschwanden (2012); Data analysis with no background subtraction are marked with the symbol *. 
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Fig. 1. — Schematic of flare parameters: flare start time t s , flare peak time t p , flare end time t e , prcflarc 
background Jbg, flux threshold f thresh, flare peak flux fp, and maximum flux-time derivative (df /dt) max . 
The curve with thick linestyle represents the rebinned and smoothed light curve which defines the local flux 
maxima and minima. 
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Fig. 2. — Example of automated flare detection of GOES observations during the day of 2012-Jan-27, which 
featured an XI. 8 flare. A number of 36 flares are detected down to the GOES B5-class level (with the starting 
times t s marked with vertical bars). Only flares with a GOES class of > CI .0 are labeled. The upper light 
curve is the 1-4 A wavelength range, which is used for flare detection, while the lower curve shows the 0.5-8 
A data. 
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Fig. 3. — Monthly flare rate during 1975-2011 on a linear scale (top panel) and on a logarithmic scale 
(bottom panel). The rate of C (blue), M (green), and X-class flares (red) are shown in colors. Note the 
proportionality of detected flares in different magnitudes. 
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Fig. 4. — Occurrence frequency distributions of GOES 1-4 A peak fluxes of solar flares detected by year. 
The powerlaw slope ap is indicated for each year. Incomplete years are the first (1974 Nov-Dec) and the 
last one (2012 Jan-Mar). 
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Fig. 5. — Occurrence frequency distribution of the GOES 1-4 A peak fluxes (top panel), the time derivative 
(df/dt) max (middle panel), and flare rise times t rise during the solar cycle maximum (years with N ev > 10000 
events), the solar cycle minimum (years with N ev < 3000 events), and for intermediate years. Note the 
invariance of the powerlaw slope for the peak fluxes (top panel) and time derivatives (middle panel) . 
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Fig. 6. — Variation of powerlaw slopes ap(i), flux-time derivative oidf/dt{t), rise time ar(t), turnover flux 
value Fi (in units of 10~ 9 W m -2 ), turnover value of time derivative (df/dt)i, (in units of 3 • 10~ n W m~ 2 
s _1 ), turnover value of rise time T\(t) (in units of s), and flare rate Nfi are (t) during the three last solar 
cycles. The theoretically predicted values are marked with thick solid lines). 
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Fig. 7. — Simulation of the fractal-diffusive SOC model for an Euclidean dimension S — 3, showing the time 
evolution of the fractal dimension Ds(t) (top panel), the instantaneous energy dissipation rate f(t) and peak 
energy dissipation ratep(t) (second panel), the total time- integrated dissipated energy e{t) (third panel), and 
the soft X-ray time profile f sxr (t) (bottom panel), which results from the convolution of the instantaneous 
energy dissipation rate f(t) (second panel) with an exponential decay function with an c- folding time of 
Tdecay = 20 (shown in insert of bottom panel) . 



